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By introducing the self-energy density functionals for the dissipative interactions between the re- 
duced system and its environment, we develop a time-dependent density-functional theory formalism 
based on an equation of motion for the Kohn-Sham reduced single-electron density matrix of the 
reduced system. Two approximate schemes are proposed for the self-energy density functionals, the 
complete second order approximation and the wide-band limit approximation. A numerical method 
based on the wide-band limit approximation is subsequently developed and implemented to simulate 
the steady and transient current through various realistic molecular devices. Simulation results are 
presented and discussed. 



I. INTRODUCTION 

Density-functional theory (DFT) has been widely used 
as a research tool in condensed matter physics, chemistry, 
materials science, and nanoscience. The Hohenberg- 
Kohn theorem [l| lays the foundation of DFT. The Kohn- 
Sham (KS) formalism Q provides a practical solution to 
calculate the ground state properties of electronic sys- 
tems. Runge and Gross extended DFT further to cal- 
culate the time-dependent properties and hence the ex- 
cited state properties of any electronic systems Q. The 
accuracy of DFT or time-dependent DFT (TDDFT) is 
determined by the exchange-correlation (XC) functional. 
If the exact XC functional were known, the KS formal- 
ism would have provided the exact ground state prop- 
erties, and the Runge-Gross extension, TDDFT, would 
have yielded the exact time-dependent and excited states 
properties. Despite their wide range of applications, DFT 
and TDDFT have been mostly limited to isolated sys- 
tems. 

Many systems of current research interest are open sys- 
tems. A molecular electronic device is one such system. 
Simulations based on DFT have been carried out on such 
devices 0, H, H, S S S M (U E3, 0- These sim- 
ulations focus on steady-state currents under bias volt- 
ages. Two types of approaches have been adopted. One is 
the Lippmann-Schwinger formalism by Lang and cowork- 
ers [7j. The other is the first-principles noneq uilibrium 
Green's function (NEGF) technique @, ©, M, El GJ, • 
In both approaches the KS Fock operator is taken as 
the effective single-electron model Hamiltonian, and the 
transmission coefficients are calculated within the nonin- 
tcracting electron model. The investigated systems are 
not in their ground states, and applying ground state 
DFT formalism for such systems is only an approxima- 
tion (lij . DFT formalisms adapted for current-carrying 
systems have also been proposed recently, such as Kosov's 
KS equations with direct current Kurth et a/.'s [l6| 
and Zheng et a/.'s [H TDDFT formulation, Cui et a/.'s 
complete second-order quantum dissipation theory (CS- 
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QDT) formalism [l8| and Burke et al.'s KS master equa- 
tion including dissipation to phonons 19]. In this paper, 
we present a new DFT formalism for open electronic sys- 
tems, and use it to simulate the steady and transient 
currents through molecular electronic devices. The first- 
principles formalism depends only on the electron density 
function of the reduced system. 

As early as in 1981, Riess and Munch [20j discovered 
the holographic electron density theorem which states 
that any nonzero volume piece of the ground state elec- 
tron density determines the electron density of a molec- 
ular system. This is based on that the electron den- 
sity functions of atomic and molecular cigenfunctions are 
real analytic away from nuclei. In 1999 Mezey extended 
the holographic electron density theorem j21j. And in 
2004 Fournais et al. proved again the real analyticity of 
the electron density functions of any atomic or molecular 
eigenstates [22| , Therefore, for a time- independent real 
physical system made of atoms and molecules, its elec- 
tron density function is real analytic (except at nuclei) 
when the system is in its ground state, any of its excited 
eigenstates, or any state which is a linear combination 
of finite number of its eigenstates; and the ground state 
electron density on any finite subsystem determines com- 
pletely the electronic properties of the entire system. 

As for time-dependent systems, the issue was less clear 
until recently we [Hf were able to establish a one-to- 
one correspondence between the electron density function 
of any finite subsystem and the external potential field 
which is real analytic in both i-space and r-space. For 
time-dependent real physical systems, we have proved the 
following theorem: [23| 

Theorem: If the electron density function of a real fi- 
nite physical system at to, p(r,io), is real analytic in r- 
space, the corresponding wave function is $(£0)7 and the 
system is subjected to a real analytic (in both i-space 
and r-space) external potential field v(r,i), the time- 
dependent electron density function on any finite sub- 
space D, /?£>(r,i), has a one-to-one correspondence with 
v(r, t) and determines uniquely all electronic properties 
of the entire time-dependent system. 

According to the Theorem, the electron density func- 
tion of any subsystem determines all the electronic prop- 
erties of the entire time-dependent physical system. This 
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FIG. 1: Schematic representation of the experimental setup 
for quantum transport through a molecular device. 



proves in principle the existence of a TDDFT formalism 
for open electronic systems. All one needs to know is the 
electron density function of the reduced system. 

This paper is organized as follows. In Sec. [IT] we de- 
scribe a TDDFT formalism for open electronic systems 
based on an equation of motion (EOM) for the reduced 
single-electron density matrix. By utilizing the holo- 
graphic electron density theorem, the self-energy func- 
tionals with explicit functional dependence on the elec- 
tron density of the reduced system are introduced, and 
thus a rigorous and efficient first-principles formalism for 
the transient dynamics of any open electronic system is 
established. Two approximate schemes, the complete 
second order (CSO) approximation for the dissipative in- 
teraction and the wide-band limit (WBL) approximation 
for the electrodes, are proposed for the self-energy func- 
tionals in Sec. |TTJ To demonstrate the applicability of 
our first-principles formalism, TDDFT calculations are 
carried out to simulate the transient and steady current 
through realistic molecular devices. The detailed numer- 
ical procedures and results are described in Sec. IIVI Dis- 
cussion and summary are given in Sec. |Vj 



II. FIRST-PRINCIPLES FORMALISM 

A. Equation of motion 

Fig. [T] depicts an open electronic system. Region D is 
the reduced system of our interests, and the electrodes 
L and R are the environment. Altogether D, L and R 
form the entire system. Taking Fig. [T] as an example, 
we develop a practical DFT formalism for the open sys- 
tems. Within the TDDFT formalism, a closed EOM has 
been derived for the reduced single-electron density ma- 
trix a(t) of the entire system [24 |: 



i&(t) = [h(t),a(t)}, 



(1) 



where h(t) is the KS Fock matrix of the entire system, 
and the square bracket on the right-hand side (RHS) de- 
notes a commutator. The matrix element of a is defined 
as &ij (t) = (a,j(t) ai(t)), where ai(t) and flj(t) are the 
annihilation and creation operators for atomic orbitals i 
and j at time t, respectively. Fourier transformed into 
frequency domain while considering linear response only, 



Eq. ([1} leads to the conventional Casida's equation [2o| . 
Expanded in the atomic orbital basis set, the matrix rep- 
resentation of a can be partitioned as 



cl &ld <Jlr 

O'DL &D &DR 
0~RL <?RD CR 



(2) 



where ox, an and o~d represent the diagonal blocks cor- 
responding to the left lead L, the right lead R and the de- 
vice region D, respectively; ctld is the off-diagonal block 
between L and D; and ord, &lr, &dl, ^dr and cjrl are 
similarly defined. The KS Fock matrix h can be parti- 
tioned in the same way with a replaced by h in Eq. ([2]). 
Thus, the EOM for ctd can be written as 



ia D = [h D ,a D ]+ (h Da a aD — a Da h aD ) 

a=LM 

= [ko, <?X>] - i Qa, 
a=L,R 



(3) 



where Ql (Qr) is the dissipation term due to L (R). 
With the reduced system D and the leads L/R spanned 
respectively by atomic orbitals {1} and single-electron 
states {fc Q }, Eq. ([3]) is equivalent to: 



io~nm ^ ^ {J^nl^lm ^nO^lrn) i ^ ^ Qa.nmi (4) 

E-D a=L,R 

^ (h n k a O'k a m — &nk a hk a m), (5) 



i 



k a Ea 



where m and n correspond to the atomic orbitals in re- 
gion D; k a corresponds to an electronic state in the elec- 
trode a (a = L or R). h n k a is the coupling matrix el- 
ement between the atomic orbital n and the electronic 
state k a . The transient current through the interfaces 
Sl or Sft (see Fig. [T]) can be evaluated as follows, 



J a {t) 



= ~// r ^ (M) 



= (h ka i<Ji ka -<Jk a ihik a 

= -^Q a ,u = -ix[Q a (t)]. 



(6) 



l£D 



Since the dissipation term Q a {t) is not known a pri- 
ori, Eq. (J3j> is not self-closed. Therefore, at this stage 
EOM ^ cannot be solved straightforwardly to obtain 
the transient dynamics of the reduced system D. 

According to the holographic electron density theorem 
of time-dependent physical systems, all physical quan- 
tities are explicit or implicit functionals of the electron 
density in the reduced system D, p (r, t). Q a of Eq. ^ 
is thus also a functional of p£>(r,i). Therefore, Eq. ([3]) 
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can be recast into a formally closed form, 



Appendix 13 cf. Eq. ©): 



l(T D 



h D [t;p D (r,t)],a D ] -i ^ Q a [t; pu(r, i)]. (7) Qc 



a=L,R 



It would thus be much more efficient integrating Eq. ([7]) 
than solving Eq. (HI), provided that Q Q [£; /Q£>(r, t)] or its 
approximation is known. We therefore have a practi- 
cal formalism for any open electronic systems. Neglect- 
ing Q a [t; p_o(r, t)] from Eq. (f7|) leads to the conventional 
TDDFT formulation [13] for the isolated reduced system, 
while Q a [t; Pd(t, t)] accounts for the dissipative interac- 
tions between D and L or R. Eq. (0) is the TDDFT EOM 
for open electronic systems, and is formally analogous to 
the master equations derived for the system reduced den- 
sity matrix in conventional QDT [26| . 

Our formalism is similar in its form to one of our 
early works, in which a dynamic mean-field theory for 
dissipative interacting many-electron systems was devel- 
oped [13, HI]. An EOM for the reduced single-electron 
density matrix was derived to simulate the excitation and 
nonradiative relaxation of a molecule embedded in a ther- 
mal bath. This is in analogy to our case although our 
environment is actually a fermion bath instead of a bo- 
son bath. More importantly, the number of electrons in 
the reduced system is conserved in Refs. [13, [H[ while in 
our case it is not. 

Burke et al. extended TDDFT to include electronic 
systems interacting with phonon baths [l9j |. they proved 
the existence of a one-to-one correspondence between 
v{v, t) and p(r, t) under the condition that the dissipative 
interactions (denoted by a superoperator C in Ref. fl9l ]) 
between electrons and phonons are fixed. In our case 
since the electrons can move in and out the reduced sys- 
tem, the number of the electrons in the reduced system 
is not conserved. In addition, the dissipative interactions 
can be determined in principle by the electron density of 
the reduced system. We do not need to stipulate that the 
dissipative interactions with the environment are fixed as 
Burke et al. And the only information we need is the 
electron density of the reduced system. In the frozen 
DFT approach [29| an additional kinetic energy func- 
tional term caused by the environment was introduced 
to account for the interaction between the system and 
the environment. This additional term is included in 
Q a [t; p D (r,t)] of Eq. ©. 



B. The dissipation term Q a 

The challenge now is to express Q a [t', Pd{y, t)]. Based 
on the Keldysh formalism |30l| and the analytical contin- 
uation rules of Langreth [3l|], Q a {t) can be calculated by 
the NEGF formulation as described in Reference [HJ (see 



leD J-oo 



G<(t,r)^ lm (r,t) + 



G r nl(t,T)X< lm (T,t)+K, 



(8) 



where G r D and are the retarded and lesser Green's 
function of the reduced system D, and and are 
the advanced and lesser self-energies due to the lead a 
(L or R), respectively. Combining Eqs. ([6|) and (jSJ), we 
obtain 



J a (t) 



2K 



dr tr 



G r D (t,T)^<(r,t) 



(9) 



Eq. ([9]) has been derived by Stefanucci and Almbladh [33| 
within the framework of TDDFT under the assumaption 
that the partitioned (34[ and partition-free [35| schemes 
are equivalent. 

It is important to emphasize that Eq. ([8]) is derived 
from the initial ground state at t — — oo when the device 
region and the leads are completely isolated, denoted by 
$o- This corresponds to the partitioned scheme devel- 
oped by Caroli et al. [34| . The dissipation term Q a can 
also be derived from the initial ground state at t = to 
when the device region and the leads are fully connected, 
denoted by ^>o, as follows (see Appendix 151 for detailed 
derivations), 

Qa, nm (t) = \c&, nm {t) / dr [g<, (t, r)S» )iro (r, t) 



leD" L o 

G r nl(t,T)E<(r,t) } + R.c, (Hi) 



where is the time immediately after to, and the first 
term on the RHS, Q Q a nm (t), arises due to the initial cou- 
plings between the reduced system and the environment. 
Eq. (flQ|) thus follows the partition-free scheme proposed 
by Cini [35[, and its associated Green's functions and 
self-energies are defined differently from those in Eq. (jSJ). 

Based on Gell-Mann and Low theorem 36], in most 
cases can be reached from <!>o by adiabatically turning 
on the couplings between the device and the leads from 
t = — oo to to. In these circumstances, the partitioned 
and partition- free schemes are formally equivalent, since 
the history of the couplings between the device and leads 
only determines V&o and its corresponding electron den- 
sity function p(r 7 to), and does not influence the dynamic 
response of the reduced system to external potentials af- 
ter to explicitly. In few cases where the turn-on of the 
couplings results in an excited eigenstate at to, Eq. ([8]) is 
only an approximation for the Q a derived from i n the 
partition-free scheme, and in principle we need to resort 
to Eq. (HUD- 
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C. Solution for steady-state current 

In cases where steady states can be reached, the 
system-bath coupling, r) = h nka (t) h kam {r) , be- 

comes asymptotically time-independent as t,r — > +00. 
The Green's functions and self-energies for the reduced 
system D rely simply on the difference of the two time- 
variables [Hj], i.e., Gd(£,t) ~ Go{t ~~ t) and E(i,r) ~ 
E(t — r), and thus we have 



E 

xE<(f 1( < 2 )G» (i 2l r) 



dta G r np {t, ti) 



* E E E 

p,q£D a=L,Rl a £a 
00 

— ieft 



X 



dfcie -te ' Sl G^(t-ti) 



OO 



- pq 



dt 2 e ie « t2 G" (t 2 -r) 



E E 



E^ 



-ief (i— r) 



p,q£D a=L,R l a £a 



(11) 



G r /(e) 

yir,a 
nm 



(0 



[e/-/i D (oo)-S'-' (e)]- 1 , (12) 

E Er-^-^^r 1 ' ( 13 ) 



a=L,R, lea 



where / is an identity matrix, 6 is an infinitesimal positive 
number, and ff is the occupation number of the single- 
electron state l a of the isolated lead a (L or R). The 
steady-state current can thus be explicitly expressed by 
combining Eqs. JTT|> - (fT3)l . 



J L (oo) = 



-Jr(oo) 

- QL.nn(oo) 



l fcGL lER 

xti-[G r D (e^r l -G a D (e^)T k 

-E/*E^- e ?) 



x tr 



R\ r l 



= J [f L (e)-f R (e)]T(e)de, 



(14) 



T(e) = 2vr ?7i ^tr[G^(£)r«(6)G? 5 (6)r i (6)j.(15) 

Here T(e) is the KS transmission coefficient, /"(e) is the 
Fermi distribution function, and rj a (e) = X^fcea — e fc) 



is the density of states (DOS) for the lead a (L or R). 
Eq. (I14p appe ars formally analogous to the Landauer for- 
mula [23, [39| adopted in the conventional DFT-NEGF 
formalism (91 ITlT] . However, to obtain the correct steady 
current, the noncquilibrium effects need to be properly 
accounted for. This may be accomplished by substitut- 
ing the asymptotic values of the TDDFT XC potential 
for the ground state DFT counterpart in Eq. (fT4")l . 



D. Self-energy functionals 

Due to its convenience for practical implementation, 
Eq. ([8]) is adopted in our formalism. The Green's func- 

ij can be calculated via the 
following EOMs if S° and E< are known, 

: dG r nm (t,r) 



tions G r D and G^ in Eq. 

E° and E< are : 

6(t - t) S nm + Ki(t) G r lm {t, t) 

dt^ nl (t,t)G r lm (i,r), (16) 



Of 



E 

leD 



r) G < 



(t,r) 



Of 



E 



dt 



E<(f,t)GL(t» + S r ni (M] 



x Gf m (t, r) + V h nl (t) G< m (t, r), (17) 



leD 



Ea=L,/?. S a, and 



where 17 = E a=L , R (K)^ ^< 
Gfj = (G r D y . The key quantities for the evaluation 
of Q a [t; Pd(y, t)} are thus the self-energies S° and E<. 
According to our Theorem, E° and E< are in princi- 
ple functionals of pn(r,t). Therefore, instead of find- 
ing Q a [t; pd(j, t)} directly, we need now seek for the 
density functionals E£[t, t; pu(r, t)] and E<[t, t; pr>(r, t)]. 
By their definitions, the self-energy terms have explicit 
functional dependence on the electron density function 
of the entire system, p = {pd, p a ): 



id(t — t) h[r; p] exp J h a [t; p a ]dt 
x h[t;p], 

ih[r;p] f a (h a [-oo; p a ]) 



K[r,f,p\ 

s a b~, t; p] 



where p a is the electron density function in the 
h a is the KS Fock matrix of the isolated lead a, 
is the Fermi distribution function for a (L or R) 
on our Theorem, p a are determined uniquely by 
a certain continuation (CT) operation, i.e., 



(18) 



xexplil h a [t; p a ]dt > h[t; p], (19) 



lead a, 
and f a 
Based 
Pd via 

(20) 
(21) 

(22) 
(23) 



CT 



Pc{r,t), 
<' T [v,t; p D {v,t)\. 



PD{r, t) 

p a (r,t) = p~ 
We obtain thus the following functionals, 

E£(r,t) = K [T,t;p D ,pZ T [p D ]] , 
E<(r,i) = E< [r,t;p D ,p^ T [p D ]] . 
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FIG. 2: The molecular device region D is subject to the 
boundary conditions AV L (t) and AV R {t) at the interfaces 
Sl and Sr. The interactions between the region D and the 
lead L and R are accounted for by the self-energy functionals 
El and E_r, respectively. 



Note that the CT operation is case dependent, and of- 
ten approximate in practice. For the system depicted in 
Fig. [TJ the CT operation from pu to p a may be approx- 
imated by a translation over repeating unit cells if the 
bulk electrodes are periodic, i.e., 



Pa{r, t) 



Pa T M «pu(r + JVR,t), 



(24) 



where t = refers to the initial time when the entire 
connected system is in its ground state, R is the base 
vector perpendicular to the interface S a for the lead a, 
and N denotes an integer which makes the translated 
vector r + AR to be inside the reduced system D as 
well as near the interfaces S a . To ensure the accuracy of 
such an approximate CT operation, it is vital to include 
enough portions of electrodes into the region D, so that 
the electron density function near the interfaces S a takes 
correctly the bulk values. 

Of course, there could be cases that the approximate 
Pu T [pd} may deviate drastically from their exact val- 
ues some distance away from the boundary. Usually 
E° and £< depend mostly on the electron density near 
the boundary where the approximate Pa T [pD\ agree best 
with the correct p a . The resulting S°[prj, Pu T [pd\] and 
[pDi Pa T [pd]] thus provide reasonable approximations 
for their exact counterparts. For cases where the self- 
energies happen to rely heavily on p a far away from D, 
the approximated CT breaks down, and our method fails 
to be applicable. 

Given T, a [pn] and £<[/?£>] how do we solve the 
EOM in practice? Again take the molecular device 
shown in Fig. [T] as an example. We focus on the re- 
duced system D as depicted in Fig. [21 and integrate the 
EOM ([7]) directly by satisfying the boundary conditions 
at Sl and Sr. AV L (t) and AV R (t) are the bias voltages 
applied on L and R, respectively, and serve as the bound- 
ary conditions at Sl and Sr, respectively. At t — > — oo, 
AV L = AV R = 0, and AV L (t) and AV R (t) are turned 
on near t = 0. We need thus integrate Eq. J7J) together 
with a Poisson equation for the Coulomb potential inside 
the device region D subject to the boundary condition 
determined by the potentials at Sl and Sr. It is impor- 
tant to point out that Q a [t; Pd(y, t)] is actually a nearly 



local quantity of the reduced system through the local 
coupling matrix terms hn a (a = L or R). In this sense, 
our formalism for open electronic systems is not in con- 
flict with the " nearsightedness" concept of Kohn [37j • 



III. TWO APPROXIMATE SCHEMES FOR 
SELF-ENERGY DENSITY FUNCTIONALS 

A. Complete second order approximation for 
dissipative functional 

Eqs. ([5]) and (flTj]) appear quite complicated. To 
have an unambiguous interpretation of the dissipation 
term Q a , we further assume the KS Fock matrix hr> is 
time- i ndep endent and treat G^(f, t) by means of CS- 
QDT 26]. Eq. © is thus simplified to be (see Ap- 
pendix [C] for details) 

Q a (t) = i{[£>(h D ),cr D y + [E<(M»^] + } , (25) 

where <jd = I — &d is the reduced single-hole density ma- 
trix of the reduced system. On the RHS of Eq. (|25p a new 
commutator has been introduced for arbitrary operators 
A and B: 



[A, B]t = AB - B^AK 



(26) 



E„ :> are the causality-transformed counterparts of ^ , 
with £<>>(£, r) = — r) presumed, i.e., 



£<'>(h D ) = / die*"* £<•>(«) 
Jo 

= TT a ±) (h D )±iA a ±) (h D ), (27) 

where 1^ (/id) and A. a (hn) are real symmetric ma- 
trices, and associated with each other via the Kramers- 
Kronig relation 26]. Therefore. Eq. (|2"5)) can be expanded 
as 

Q a (t) = i [r(-\h D ),a D ~\ + {A a -\h D ),a D } - 

r£\h D ),v D ]-{AM{h D ),it D }. (28) 



The physical meaning of Eq. (|28|) is clear and intuitive: 
the first and third terms on its RHS account for the en- 
ergy shifting of occupied and virtual orbitals of the re- 
duced system due to the couplings with the lead a, re- 
spectively; and the second and fourth terms on its RHS 
are responsible for the level broadening of occupied and 
virtual orbitals in D due to the lead a while contributing 
to the transient current, respectively. The second term 
accounts for the electrons leaving the device region, and 
the third term describes that the holes hop onto the elec- 
trodes or the electrons enter the device region from the 
electrodes. 



B. Solution for transient current with WBL 
approximation and test on a model system 



JlO) 



J R (t) 



To simplify the solutions of Eqs. ([Io ) -([rT l) . the WBL 
approximation (32l . |40T | may be adopted besides the ap- 
proximate CT operation (c/. Eq. (|24|l ). which involves 
the following assumptions for the leads: (i) their band- 
widths are assumed to be infinitely large, such that 
the summation over all the single-electron states in the 
leads can be replaced by an integration over the en- 
tire energy range, i.e., J^kca ~* f-oo der )c*( e )> (h) their 
line- widths, A^(i, t), defined by the DOS at Sl or 
Sr times the system-bath couplings, i.e., A^(t, r) = 
7T77 Q (e^) T kc ' (t, t), are treated as energy independent, 
i.e., A%(t,r) « A Q (t,r) w A Q , and (iii) the level shifts 
of L or R are taken as a constant for all energy levels, 
i.e., Aejj*(i) « Ae a (<) = -Ay a (<), where Al/ a (t) are the 
bias voltages applied on L or R at time i. 

Within the WBL approximation, the self-energy func- 
tionals can be expressed by [231 ] 



a,nm 



(r,t) = i8(t-T)A% m [p D ], 
(r,t) 



2i 


cxp /i 






7T 






■ <. + oo 


X 






.J — oo 



ic(t-T) 



(29) 



(30) 



Here AT/"(f) is the bias voltage applied on the lead a, 
and A a [p£>] is the line- width matrix due to lead a [23j |. 

KJiPd] = irri a {e f )(h nkf [p D ,p D {r + NR)} 

x h kfm [ PD ,p D (r + NK)}) , (31) 

where rj a (ef) is the density of states for a at its Fermi en- 
ergy ef, kf is a surface state of a at e/, and (• • • ) denotes 
the average over all surface states at e/. Eqs. (f29 |) — (|3T|) 
provide thus the explicit dependence of and on 
Po(r, t). 

Note that S^p^] depends on the applied voltage 
AV a (t) explicitly. In principle AV a (t) is a functional 
of po{Y,t) as well. pi)(r,i) is unknown and needs to be 
solved. The potential v(r) in DFT formalism, which in- 
cludes the potentials from nuclei and external sources, 
is a functional of electron density p(r). In any practical 
implementation of DFT, v(r) is given and used to solve 
for p(r), instead of determining v(r) from p(r). In our 
formalism AV a (t) is given as a known function and used 
to determine pi)(r,t) in the same fashion. 

Based on Eqs. (|29[) — (13~T1) . the dissipation term within 
the WBL approximation, Q^ BL , can be obtained readily 
as follows (see Appendix |D] for detailed derivations) , 



Qr L (t)=K a (t) + {A a [ PD ],a D } 



(32) 



E 
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FIG. 3: Model system for the test of the WBL self-energy 
functionals where a single site spans the device region D. 
Transient currents through leads L and R, Jh(t) and Ji?(t), 
are simulated. The inset shows the time-dependent level shift 
of lead R. 



mutator, and K a (t) is a Hermitian matrix, 



+ 1 [i-u a (t) 

de 



e iet l x 



— - A« + H.c,(38) 



Here the curly bracket on the RHS denotes an anticom- 



e-h D (t) + iA + Ae a (t) 



where pP is the chemical potential of the entire system, 
the overall line- width A = ^ Q A", and the effective prop- 
agator of the reduced system U a (t) is 

U a (t) = e _i /o t [' iD ( r )- iA - AeQ < v )] dT . (34) 



C. Numerical test of wide- band limit 
approximation 

The WBL approximated self-energy functionals are 
then tested by calculations on a model system which 
has previously been investigated by Maciejko, Wang and 
Guo In this model system the device region D 

consists of a single site spanned by only one atomic or- 
bital (see Fig. [3]). Exact transient current driven by a 
step voltage pulse has been obtained from NEGF simu- 
lations [40], and the authors concluded that the WBL 
approximation yields reasonable results provided that 
the band-widths of the leads are five times or larger 
than the coupling strength between D and L or R. 
The computational details are as follows. The entire 
system (L + R + D) is initially in its ground state 
with the chemical potential pP. External bias voltages 
are switched on from the time t = 0, which results 
in transient current flows through the leads L and R. 
Sh D (t) =h D (t) -h D (0), Ae L (t) and Ae R (t) are the level 
shifts of D, L and R at time t, respectively. In our works 
we take Sh D (t) = \ [Ae L (t) + Ae R (t)] , Ae L (t) = 0, and 
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FIG. 4: The calculated transient current through Sr within 
the WBL scheme. We set fx — hn(0) = for the ground 
state; and Ae L (t) = and Ae R (t) = Ae R (1 - e~ t/a ) after 
switch-on. The above panels show different cases where (a) 
Ae R = 2 eV, A L = A R = 0.1 eV; (b) Ae R = 0.2 eV, A L = 
A R = 0.1 eV; (c) Ae R = 10 eV, A L — A R = 0.1 eV; and (d) 



2 eV, A 



A 



0.04 eV, respectively. 



Ae R (t) 



Ae R {l 



-t/a 



), where a is a positive con- 



stant. The real analytic level shift Ae R (t) resembles per- 
fectly a step pulse as a — > + (see the inset of Fig. [3]). 
The calculation results are demonstrated in Fig. 2J We 
choose exactly the same parameter set as that adopted 
for Fig. 2 in Ref. [4(|, and the resulting transient cur- 
rent, represented by Fig. [Ha), excellently reproduces the 
WBL result in Ref. [40(, although the numerical proce- 
dures employed are distinctively different. The compar- 
ison confirms evidently the accuracy of our formalism. 
From Fig. BJa) — (c) it is observed that with the same 
line- widths A Q , a larger level shift Ae^ results in a more 
fluctuating current, whereas by comparing (a) and (d) we 
see that under the same Ae R , the current decays more 
rapidly to the steady state value with the larger A a . 

By transforming its integrand into a diagonal repre- 
sentation, the integration over energy in Eq. (|33p can 
be carried out readily. Therefore, Q^f BL are evaluated 
straightforwardly, which makes the above solution pro- 
cedures for transient dynamics within the WBL approx- 
imation a practical routine for subsequent TDDFT cal- 
culations. 



IV. TDDFT CALCULATIONS OF TRANSIENT 
CURRENT THROUGH MOLECULAR DEVICES 

A. Numerical procedures 

With the EOM and the WBL approximation for 
the self-energy functionals and £<[/5d], it is now 

straightforward to investigate the transient dynamics of 
open electronic systems. All our first-principles calcu- 



lations are carried out with a self-developed package 
LODESTAR [HI]. 

The ground state properties of the reduced system 
at t — are determined by following the partitioned 
scheme approach adopted in conventional DFT-NEGF 
method [g, [^, [l(| EH- Different from the popu- 
lar periodic-boundary-condition-based approach [id . Il3l . 
HH, what we employ is a molecular-cluster-based tech- 
nique [4l[. The ground state KS Fock matrix of an 
extended cluster, covering not only the device region 
D but also portions of leads L and R, is calculated 
self-consistently by conventional DFT method with local 
density approximation (LDA) for the XC functional [2]. 
Its diagonal blocks corresponding to the leads L and 
R are then extracted and utilized to evaluate the sur- 
face Green's function of isolated lead a (L or R), g r a = 
Q T r. \p-°\ P~ T [pp]]i by applying the translational invari- 
ance (43| (c/. Eq. (f2"4"| ). In this way the possible mis- 
alignment for the chemical potentials of the isolated leads 
L and R, especially when they are made of different ma- 
terials, can be avoided so long as the extended cluster 
is chosen large enough. In an orthogonal atomic orbital 
basis set, the line-widths A a [p£>] within the WBL ap- 
proximation are obtained from g r a via 



k a [ PD ] = -%{h Da g r a K;p£ T M h aD ) 



(35) 



At t = the left-hand side (LHS) of the Eq. vanishes. 
The EOM ([7]) reduces thus to a nonlinear equation for 
ctd(0), and can be solved readily by employing the NEGF 
approach as follows, 



<td(0) 



a,Q, 



(36) 



where 



G r 6°(e) = [GS°(e)]t = [e - h D (0) + iA}' 1 . (37) 



Eq. (|36|) provides the initial condition for the EOM (J7J. 
The molecular device is switched on by a step-like volt- 



age AV R (t) = -Ae R (t) = AV R (1 - e 



-t/a) 



applied on 



the right lead with a — > + (see the inset of Fig. [3]), 
while AV L (t) = 0. The self-energy functionals X£[ad] 
and E<[jOd] can be evaluated through Eqs. (f2"9")l — ([37))) 
and (|35p . The dynamic response of the reduced sys- 
tem is obtained by solving the EOM ([7]) in time domain 
within the adiabatic LDA (ALDA) Hf| for the XC func- 
tional. The induced KS Fock matrix of the reduced sys- 
tem, 5ho{t) = ho{t) — hu(0), is comprised of Hartree 
and XC components [HI , i-e., 



where 



5h D (t) = 5h%(t) + 5h% c (t), 



6hm= / dr#(r)<fo"(r,i)^(r). 



(38) 



(39) 



Here the Hartree potential 8v H (r,t) satisfies the follow- 
ing Poisson equation for the device region D subject to 
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FIG. 5: A graphene-alkene-graphene system adopted in 
TDDFT calculations. 

boundary conditions AV a (t) at every time t: 

\7 2 5v H (r,t) = -ATrSp D (r,t) 
5v H (r,t)\ SL = AV L (t) (40) 
Sv H (r,i)\ SR = AV R (t). 

To save computational resources we calculate 5hp C (t) to 
its first-order change due to the switch-on potential: 

= E V $£n knW-^(O)], (41) 
7nn£D 

x #(r)fc(r), (42) 

where u [r, i; p/j] is the XC potential. The reduced 
system is propagated from t = following the EOM Q 
by the fourth-order Runge-Kutta algorithm [44| with the 
time step 0.02 fs. Virtually the same results are yielded 
by adopting a much smaller time step, which justifies the 
accuracy of our time evolution scheme. 

B. Calculation on a graphene-alkene-graphene 
system 

A realistic molecular device depicted in Fig. [5] is taken 
as the open system under investigation. The device re- 
gion D containing 24 carbon and 12 hydrogen atoms is 
spanned by the 6-31 Gaussian basis set, i.e., altogether 
240 basis functions for the reduced system. The leads 
are quasi-one-dimensional graphene ribbons with dan- 
gling bonds saturated by hydrogen atoms, and the entire 
system is on a same plane. The extended cluster contains 
totally 134 atoms. 

In Fig. [5] we plot the calculated transient currents 
through the interfaces Sl and Sr, Jh{t) and Jn(t), un- 
der various turn-on voltages. As depicted in Fig. [13 J tit) 
and Jn(t) increase rapidly during the first few fs and 
then approach gradually towards their steady state val- 
ues. This agrees with previous investigations on model 
systems [H, |4(J. The steady currents through Sl and 
Sr are (a) -5.9 /xA and 5.9 [iA, (b) -14.2 /iA and 
14.2 /iA, (c) -18.0 ^A and 18.0 /iA, and (d) -21.3 /iA 
and 21.3 /iA, respectively, and thus cancel each other 
out exactly, as they should. By comparison of panels 




5 10 15 20 25 5 10 16 20 26 

time (Is) time {fs> 



FIG. 6: The solid (dashed) curve represents the transient 
current through the interface Sr (Sl) of the graphene-alkane- 
graphene system driven by a step-like voltage applied on the 
lead R with the amplitude (a) AV R = -0.1 V, (b) AV R = 
-0.3 V, (c) AV R = -0.5 V, and (d) AV R = -1.0 V. 



(a) — (d) it is obvious that a larger turn-on voltage re- 
sults in a more conspicuous overshooting for the tran- 
sient current. Complex fluctuations are also observed 
for the time-dependent currents, which are due to the 
various eigenvalues possessed by the nonnegative definite 
line-widths A" with their magnitudes ranging from to 
4.1 eV, corresponding to various dissipative channels be- 
tween D and L or R. From Fig. O the characteristic 
switch-on time for the graphene-alkene-graphene system 
is estimated as about 10 ~ 15 fs for applied bias voltages 
ranging from 0.1 V to 1.0 V. For much higher turn-on 
voltages the linearized form of 8h^ c (t) (Eq. (|4"Tj) ) be- 
comes inadequate, which makes such a TDDFT calcula- 
tion computationally demanding with our present coding. 

It is noted that the reduced system remains in its 
ground state in absence of an applied bias voltage. This 
is confirmed by a free propagation for the reduced sys- 
tem. During the course the transient current J/,(t) or 
Jfj(i) vanishes correctly at every time t > 0. This thus 
validates that the WBL approximated self-energy func- 
tionals derived from the partitioned scheme (c/. Eq. ©) 
is well adapted to a TDDFT formalism. 



C. Calculation on a CNT-alkene-CNT system 

The second molecular device we calculate is sketched in 
Fig. where a linear alkene is connected to semi-infinite 
single-walled carbon nanotubes (CNT) (5, 5) at its both 
ends. The device region D consists of 88 carbon and 
22 hydrogen atoms, i.e., altogether 836 basis functions 
for the reduced system. The extended cluster for the 
ground state calculation contains totally 290 atoms. The 
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FIG. 7: A CNT-alkene-CNT system adopted in TDDFT cal- 
culations. 
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FIG. 8: The solid (dashed) curve represents the transient 
current through the interface Sr (Sl) of the CNT-alkene- 
CNT system driven by a step-like voltage applied on the lead 



R with the amplitude (a) AV R = -0.1 V, (b) AV 1 



-0.3 



V, (c) AV 1 



-0.5 V, and (d) AV 1 



-1.0 V. 



calculated transient currents driven by step-like turn-on 
voltages AV R (t) (see the inset of Fig. [3]) are plotted in 
Fig. [H] Here we have set AV L = 0. The switch-on time 
for the CNT-alkene-CNT system is about 10 fs for ap- 
plied voltages ranging from 0.1 V to TO V. 



D. Calculation on an A1-C7-A1 system 

Another open system adopted in our first-principles 
calculations is depicted in Fig. [21 where a linear chain 
of seven carbon atoms is embedded between two semi- 
infinite Al leads in the (001) direction of bulk Al. The 
current-voltage characteristics of this A1-C7-A1 system 
with the same geometric configuration has been investi- 
gated extensively fiol Il3l|. In our calculation, the device 



region D consists of 7 carbon and 18 Al atoms, i.e., al- 
together 297 basis functions for the reduced system, and 
the extended cluster for ground state calculation contains 
totally 115 atoms. 

The calculated non-WBL transmission coefficient, 
T(e; AV R = 0V), is plotted in Fig. [TO] The main fea- 
tures of our result agree reasonably with those exhib- 
ited in literature [id 1 1 31 ] . The quantitative discrepancies 
may be due to the different techniques employed. For 
instance, a finite molecular cluster is explicitly treated 
in our calculation, whereas an infinite periodic system is 
considered in Refs [l(| EH, and also the basis set and 
XC functional adopted are distinctively different. The 
calculated transient currents driven by step-like turn-on 
voltages AV R (t) (see the inset of Fig. [3J are plotted in 
Fig. [TTJ The switch-on time for the A1-C7-A1 system is 
about 3 ~ 5 fs for applied voltages ranging from 0.1 V to 
0.5 V. 



V. DISCUSSION AND SUMMARY 

Kurth et al. have proposed a practical TDDFT ap- 
proach combined with the partition- free scheme (l6| . A 
number of relevant technical issues have been addressed, 
for instance, how the intractable propagation of the KS 
orbitals of an infinitely large system is transformed into 
the time evolution of KS orbitals in a finite open sys- 
tem subject to correct boundary conditions, how the 
time-dependent KS equation for the entire system is dis- 
cretized in both r and t spaces, etc.. The performance 
of their approach has been illustrated by calculations for 
one-dimensional model systems. Our first-principles for- 
malism for open electronic systems is fundamentally dif- 
ferent: (i) In our method the KS reduced single-electron 
density matrix is used as the basic variable while in 
Ref. [16( the occupied KS single-electron orbitals are 
propagated, (ii) The concept of self-energy functional is 
introduced in our formalism. In principle the self-energy 
functional depends only on the electron density function 
of the reduced system, and hence we need only focus on 
the reduced system of interest without treating explicitly 
the environment. The influence of the environment enters 
via boundary conditions and the self-energy functionals. 
This is not only for quantum transport phenomena, but 
also for any dynamic process in any open electronic sys- 
tem. In this sense we expect the EOM ([7]) to be a gen- 
eral recipe for open system problems, (iii) Our EOM 
is formally analogous to the master equations derived 
from the conventional QDT [26|. From this perspective, 
well-established methods and techniques of QDT may be 
employed to improve the evaluations of self-energy func- 
tionals and the dissipation term Q a [t; pu(r, t)] system- 
atically. For instance, another EOM has recently been 
proposed by Cui et al. based on the CS-QDT with a 
self-consistent Born approximation (SCBA) [18]. 

In conventional QDT 26] the key quantity is the re- 
duced system density matrix, whereas in Eq. ([7]) the ba- 



10 




FIG. 9: A linear carbon chain is sandwiched between two Al 
leads in the (001) direction of bulk Al. 
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of electron density only, Exc = Exc[p{ r )]- This is true 
when the interaction between the electric current and 
magnetic field is negligible. However, in the presence 
of a strong magnetic field, Exc = Exc[p{ r ) Jp( r )] or 
Exc — Exc[p( r ), B(r)], where j p (r) is the paramag- 
netic current density and B(r) is the magnetic field [47j |. 
In such a case, our first-principles formalism needs to be 
generalized to include j p (r) or B(r). Of course, j p (r) or 
B(r) should be an analytical function in space. 

To summarize, we have proven the existence of a first- 
principles method for time-dependent open electronic 
systems, and developed a formally closed TDDFT for- 
malism by introducing the concept of self-energy func- 
tionals. In principle the self-energy functionals depend 
only on the electron density function of the reduced sys- 
tem. With an efficient WBL approximation for self- 
energy functionals, we have applied the first-principles 
formalism to carry out TDDFT calculations for transient 
current through realistic molecular devices. This work 
greatly extends the realm of density-functional theory. 



FIG. 10: Non-WBL transmission coefficient T(e; AV R 
of the A1-C7-A1 system. 
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APPENDIX A: DERIVATION OF EQ. © WITH 
THE KELDYSH FORMALISM 
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FIG. 11: The solid (dashed) curve represents the transient 
current through the interface Sr (Sl) of the A1-C7-A1 system 
driven by a step-like voltage applied on the lead R with the 
amplitude (a) AV R = -0.1 V, and (b) AV R = -0.5 V. 



sic variable is the reduced single-electron density ma- 
trix, which leads to the drastic reduction of the de- 
grees of freedom in numerical simulation. Linear-scaling 
methods such as the localized-density-matrix (LDM) 
method 0, EH may thus be adopted to further speed 
up the solution process of Eq. 0. Therefore, Eq. (0 
provides an accurate and convenient formalism to inves- 
tigate the dynamic properties of open systems. 

It is worth mentioning that our first-principles method 
for open systems applies to the same phenomena, prop- 
erties or systems as those intended by Hohenberg and 
Kohn [l|, Kohn and Sham Q, and Runge and Gross 
i.e., where the exchange-correlation energy is a functional 
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In the Keldysh formalism [30|, 
single-electron Green's function G ka 



the nonequilibrium 
m {t,t') is defined by 



G Km (t,t') = -i{T c {a ka (t) al(t')}) 



(Al) 



where Tc is the contour-ordering operator along the 
Keldysh contour [3(J [HJ (see Fig. [T2|) . Its lesser com- 



ponent, G% m (t,t'), is defined by 

l (t,t')=i(al(t')a ka (t)). 



(A2) 



The formal NEGF theory has exactly the same struc- 
ture as that of the time-ordered Green's function at zero 
temperature [32l . [48j. Thus, the Dyson equation for 
Gk a m{t,t / ) can be written as 

Gk am (t,t') = / dr g ka (t,r) h ka i(T)Gi m (T 7 t'), 
ieD ^ c 

(A3) 

where Gi m (r,t') and g ka (t,r) are the contour-ordered 
Green's functions for the reduced system D and the iso- 
lated semi- infinite lead a (L or i?), respectively, and the 
integration over r on the RHS is performed along the 
entire Keldysh contour (see Fig. [P2)l . 
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FIG. 12: The Keldysh time contour on which nonequilibrium 
Green's function theory is constructed. On the contour, the 
time T\ is earlier than ti even though its real-time projection 
appears larger. 



G r lm {r,t% Gf m (r,t') and G< n ( T ,t') denote the re- 
tarded, advanced and lesser components of G/ m (r, £'), 
respectively. Their definitions are as follows, 

Gf m (r,i') = -^(r-f'XHr),^')}), (A4) 

Gt m (r,t') = i#(t'-T)({ai(T),al(t')}), (A5) 

G< m (r,t') ee i^COoiW). (A6) 

where $(t — £') is the Heaviside step function, and the 
expectation values (• ■ • } are taken at the ground state 
of the entire system at t = — oo, i.e., when the reduced 
system and the environment are completely decoupled. 
G[ m (r, t') and Gf m (T,t') are to be calculated via their 
EOMs 0^ -CCf]). The related self-energies ££(t,r) and 
r) are evaluated through 



E &j*-(*)sfc a Mfck.»M. ( A7 ) 

E hik a (t)g<Jt,T)h kan (T), (A8) 



for a = L or i?. Here g k (t, r) and g^p (t, r) are the ad- 
vanced and lesser surface Green's functions for the iso- 
lated lead a (L or R) [32[. 

Applying the analytical continuation rules of Langreth 
[3ll |. we have 

G< k Jt',t) EE 

= -[G£ m (t,0]* 

/OO 
*"Wr) ^(r,t)G^(i',r) 
lfcJJ -°° 

+ gUr,t)G< ll {t',T)] (A9) 
by adopting the following equalities: 

G r ml {t',T) = [G? m (T,t')]* , 

G<tf,r) = -[G< n (r,t')Y, 

9l(r,t) = [g r ka {t,r)]* , 

9t(r,t) = -[g&(t,T)]\ (A10) 

Note that a mka (t) is precisely the lesser Green's function 
of identical time variables, i.e., 

^*«,(*) = -*G< jko (*,t / )| t , =t . (AH) 

By inserting Eqs. jA9]) and (|A11|I into Eq. ©, Eq. © 
can be recovered straightforwardly. 



APPENDIX B: THE DISSIPATION TERM Q Q IN 
PARTITION-FREE AND PARTITIONED 
SCHEMES 

For brevity, £ Q=Lifl £ fc g Q wil1 be shortened to J2 ka - 
The Hamiltonian of the entire noninteracting KS system 
is 

H(t) = E h rnn(t)al n a n +'^2e ka (t)al a a ka 

E E [Kike (M n ak a + H.c] . (Bl) 



mG-D fe a 

Initially (at t = t a ) the entire KS system is in its ground 
state V(t ) (denoted by |0) hereafter), i.e., H(t Q )^(t Q ) = 
-E'O^'(io)- We define the following Heisenberg creation 
and annihilation operators (h = 1): 

t fj.\ i It H(t)cIt f -iff H(t)cIt 

<W = e Jt o <e J *o 

» f/ HMdr -i A* H(r)dT 

a m {t) ee e Jt o w a ro e J *o , 
4 a (t) ee e'^^e"^"^, 
o fca (t) ee e^^^e-*^^*, (B2) 
which satisfy their respective EOMs (dt = §{)'■ 

dtal(t) = i^a\{t)h irn {t)+i^al a {t)h kaTn {t), 

ieD k a 

d t a m (t) = h mi (t)a,i(t) - i^h mka (t)a ka (t), 

d t a ka (t) = -i^2h ka i(t)a,i(t) - i e ka (t)a ka (t), (B3) 

with the initial conditions: a^io) = a mi a m(^o) = a m, 
a iL(*o) = a l a , and a ka (t a ) = a ka . 

The retarded, advanced and lesser surface Green's 
functions for the isolated lead a (L or R) are defined 
as follows, 

g%(t,r) ee T i#(±t T T)(a\{b k Jt),bijT)}\a),m 
g<(t,r) ee i(a\bUr)b ka (t)\a), (B5) 



where the curly bracket on the RHS of Eq. (|B4|) denotes 
an anticommutator, and \a) is the ground state wave- 
function corresponding to the initial lead Hamiltonian 
H a (to): 



a ka 



(B6) 



k£a 



The Heisenberg operators in Eqs. (|B4|1 and (|B5|) are de- 
fined by 

,+ /,n if* H a ( T )dT f -if? H a {r)dT 

b k (t) ee e ° aj, e Jt o , 

i /,n if' H a (r)dT -if* H a (r)dr /T} -^ 

b ka {t) ee e ° afc a e Jt ° . (B7) 
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We then define the retarded, advanced and lesser 
Green's functions of the entire KS system via their ma- 
trix elements as follows, 

G r t f(t,r) ee T ^(±tTT)(0|{a i (t),at(r)}|0>,(B8) 

G< a] {t,r) ee i(0\a](r)a ka (t)\0), (B9) 

G<^(t,r) ee i(O\al (r)a ka (t)\O), (BIO) 

G<(t,r) ee i(0\a](r)a z (t)\0). (Bll) 

where /3 = L or i?, and p,g denotes a single-electron state 
in the lead (3. Hereafter we only solve the Green's func- 
tions for time variables t and r ranging from ig to +00. 
Taking the first-order time derivatives of g k (t, r) and 
9k a (*> T ) leads to 

I^-e*.(t)]5lk.(^) = (B12) 
-[^ T + e fea (r)] 5 L(t,r) = <5(t-r), (B13) 



with the initial conditions for Eg. (|B12|) : g£ (t,r)| t=T + = 
~h 9l a {t, r)\t=r = and gl a (t,T)\ t=T - = 0; and for 
Eq. (Uni: g% a {t,r)\ T = t - = 0, ffg o (t,T)| T=t = |, and 
fffe a (*,T")|T=t+ = respectively. g r ka (t,r) and g%Jt,r) 
can thus be utilized to solve partial differential and 
integro-differential equations. For instance, we have the 
EOM for G% aj (t,T) as follows, 

[id t -e ka {t)]G< aj {t,r)= h kam {t)G< 3 {t,r). 

meD 

(B14) 

Combining Eqs. (|BT4|) and (|BT2"|) . we obtain 

G L-(*' T ) = E f + dtgl a {t,T)h kam {T)G< 3 {lr) 



+ igUt,t Q )G<(t+,r). 



(B15) 



With a similar but slightly more tedious treatment for 
the time variable r, we arrive at 



~i E G tm(t,t+)G a mj (t ,T) 
meD 

x g» (to, t) ^ (t) G" (t, r), (B16) 



where is short for X)/3=£ i? Tlpep ■ taking t = 
Eq. (|B16|) and then insert it into Eq. (jB 15|) . we have 



meD 5 o 



+ i 



E 9L( i .*o)a fcQm (t+)G I %.(to,r) 



m£D 



+ «EE /+ rf *fL(*'*0) <J fe«P/3(*0 ) 
P/3 meD l Q 

x < a (*o, t) h PBm {t) G a rn At, r), (B17) 



where the following equalities have been adopted: 

o* a m(t) = -*G< m (t,r) , (B18) 

T — t 

vk aP0 (t) = -iGl P0 (t,r)\ T=t . (B19) 
From Eq. ([5]) the dissipative term Q a (i) is expressed by 
Qa,ij (t) = h ika (t) a kaJ (t) + H.c. 

kea 

= E h ^ (*) G iw (*> T ) + H - c - ( B2 °) 

fee a 

Combining Eqs. (jB17|) and JB20}, we have thus 



Ql,H (*) + E / + d * S a,>m(*. *) *) 

jT*^ iiro (t,t)G^(t,t)| 



meD 
H.c 



(B21) 



where 



->< 

J a,im 



a,ij W 
(*,*) 



*E E ^(*)ffL(*'*o)<^ a m(*o) 

xG a mj (t Q ,t), (B22) 

« E E ^ fe ° (*) L (*' *o) CT *cPi3 (*o ) 

tea P/3 

xg; (t o ,t)h ppm (t), (B23) 
£-(*.*)'**-"»(*). (B24) 

tea 



where and Y7 a are the lesser and retarded self-energies 
of the device region, respectively. Note that by defi- 



nitions CT D = [G a D }\ G< = -[G<]+, = [E£]t and 
£< = -[S<]t, therefore, it's trivial to validate Eq. (|B2"T|) 
is equivalent to Eq. (jTTJjl . 

It is important to emphasize that ^(tj) may be differ- 
ent from ^ (to); so that the corresponding reduced single- 
electron density matrix cr(t^) may also differ from cr(to). 

This would happen if J t *° H(t) dr ^ 0, for instance, in 
the cases where the external field involves a Delta func- 
tion switched on at to. However, for real physical sys- 
tems, the applied external field is real analytic in time. 

In this circumstance, J t ° H(t)cLt = 0, *(to) = *(*o), 

and c(io") = (J {to). 

The above derivations follow rigorously the partition- 
free scheme, since the initial state "J (to) can be the 
ground state of the fully connected entire system includ- 
ing the device region and the leads. As for the partitioned 
scheme, we need to introduce another reference state $0, 
which is the ground state of Hamiltonian H, 



H = h mn {t )a m a n + ^ e ka (t )a' k 



k a a k a 



(B25) 
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Since H does not contain any coupling terms between 
D and L or R, <J> depicts the scenario that the device 
region and the leads are isolated from each other. Hence 
there is no electron populated across the boundary Sl 
and Sr. i.e., <7Da — and <tlr = 0. We now assume 
if? (to) can be reached by a time propagation of the entire 
system starting from the state $o, i.e., 



(B26) 



At t — —oo, H(—oo) = H and <r(— oo) = a. In this 
sense, the initial time for the Heisenberg creation and 
annihilation operators defined in Eq. (IB2j) becomes — oo 
instead of to, and the above derivations for the various 
Green's functions remain valid. Note that for the decou- 
pled ground state <E>o, we have 

a 4 ,(-oo) = a%, (B27) 
a kaj (-oo) = 0, (B28) 
o-fc QP(3 (-oo) = 6 a p6 kp f° a , (B29) 

where f9 is the initial occupation number of the single- 
electron state k a . Thus the Green's functions and self- 
energies previously derived can be simplified as follows, 



>(*,*) = 



G<(t,r) 



: E /*- W C (*. *0) (*0, t) ^m(t), 

E ^fc- (*) dt (*> *) ^m(t), (B30) 



meD 



+ 



dtg r ka (t, t) h kam (i) G<Ai, t) , (B31) 



' E 

rnn^D 



Gl m (t,-oc)a° mn G a nj (-oo,T) 



+ E 

rnn^D 

x E!S 



r/fi 



d*2 GTm(*)*l) 



^i^G^r). 
The dissipative term Q Q is thus expressed as 



(B32) 



E 



dt 



t (M~)G" (M) 



+ SE im (t,t)G<,(t,t)+H.c. . (B33) 



Eqs. ([B30|> — (fH33j) recover exactly Eq. flgj) derived from 
the Keldysh NEGF formalism [30(. Therefore, we con- 
clude that as long as the relation (|B26|) holds, the 
partition- free and the partitioned schemes of NEGF yield 
exactly the same dissipation term Q a (t) for t ^ to- 
la fact, Eq. (526)1 can be proved by Gell-Mann and 
Low theorem (1951) (36|, which basically states that 
if?(to) can be reached from $o by adiabatically turning on 
the coupling terms between D and L or R from t = — oo 
to to. The resulting if? (to) is an eigenstate of the Hamil- 
tonian H(to) and in most cases is the ground state. 



APPENDIX C: DERIVATION OF EQ. ([25]) 



The greater Green's function for the reduced system, 



Gij(t,r), is defined as 



G>(t,T)^-i(a % (t)a](T)). 



(CI) 



The advanced Green's function of the reduced system can 
thus be expressed as 

G a D (t,r) = -tf(r-i) [G>(t,r) - G<(t,r)] . (C2) 

Similarly the retarded and advanced self-energies can be 
associated with the greater and lesser self-energies as fol- 
lows, 

E5°(t, r) = ± 0(± * T r) [E> (t, r) - £< (f , r)] , (C3) 
where the greater self-energy £>(i, t) is defined as 



^k a (t)9 k (t,T)h ka3 (T) 



= -i} ,h ika (t)h ka j(T) 

k£a 

x(a\b ka (r)bl(t)\a). (C4) 
Eq. ([5]) is thus equivalent to 

Q a (t) = f dr[E>(t,r)G<(r,t) 

- E<(t,r)G>(r,t)] +H.c. (C5) 

In cases where the KS Fock matrix of the reduced system, 
ho, is time- independent, the greater and lesser Green's 
functions can be approximated by QDT perturbatively 
to complete second order (26|, i.e., 

G>(r,t) « e ihD ^G>(t,t) = H)^-^) a D , 
G£(r,t) « e l/lD(t ~ T) G£(M) ^je'' 10 ^^, (C6) 

where <td = I — ajj is the reduced single-hole density 
matrix of the reduced system. Assuming the lead Hamil- 
tonian to be time- independent, we have £^ :> (t, r) = 
— t). Hence, Eq. (|C5[) can be recast into 

roc 

Q a (t) = / dTE<(T)e l,lDr G>(t,i) 
Jo 

poo 

+ / drE>(r)e' i ' iDr G<(t,t)+H.c..(C7) 
Jo 



To evaluate the causality transforms involved in Eq. (|C7[) . 
we define 

1 Z" 00 

^'(M = ± T{ / dr[E<'>(T)e^ T 
z * Jo 

+ e- ih ^<'>(-r)], (C8) 
TtHhn) = T2 ^ dr[E<^(r)e^ T 

_ e -ih fl r E <,>(_ r j] _ (C9) 
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Here the equality [S^' > (t)]t = — £< :> (— r) has been 
adopted. With Eqs. ([€? !) -(|Cl) ]) . Eqs. p5 ]) -([25 1) are read- 



ily recovered. Generally A^^/i^) and T^^d) are Her- 
mitian matrices, and associated with each other via the 
Kramers-Kronig relation [26j . In particular, when the 

KS Fock matrix h is real, A^^/id) and r„ (/id) be- 
come real symmetric matrices. With Q a expressed by 
Eq. P5|) . the EOM for <td is reformulated as 



a=L,R 

+ ]T [E<(^)^D] f . (CIO) 



a=L,R 



Eq. (jCTO)) resembles closely Eq. (8) in Ref. [18|, which 
is developed from CS-QDT with the Markovian approx- 
imation. The correlation functions of the leads used 
m Ref. [3, C { a ±} (t,T), are related to the self-energies 
adopted in our work as follows, 

X<<>(t,T)=±iCtHt,r). (Cll) 

Following the SCB A scheme proposed in Ref. [18| , higher 
order effects due to interactions between the reduced sys- 
tem and the environment can be partially accounted for 
by substituting in Eq. (|C6|) an effective propagator of the 
reduced system, & lh o C*- 1 "), for the propagator of the iso- 
lated reduced system, e lh o(t- T ) ^ w here h e ^ is some effec- 
tive KS Fock matrix of the reduced system. This results 
in self-energy terms T,^^ (h^) instead of in 
Eq. dCTOl) . 



APPENDIX D: WIDE-BAND LIMIT SCHEME 
FOR THE DISSIPATION TERM Q Q 

With the WBL approximation, the advanced self- 
energy becomes local in time [3^ . 



fco£a 



hnk a (r) hk a m (t) gl ( T , t) 



^ h nka {r) h kaTn (t) 

[M(t-T)e i€ ^- T )e^ Aea © d * 



= -&(t 



e Ht~r) de I A c 



i5{t - t)K, 



(Dl) 



Here the Dirac Delta function on the RHS effectively re- 
moves the tricky off-diagonal elements of (t, r) from 
the NEGF formulation for Q a (cf. Eq. ©). The third 
equality of Eq. (|D1[) involves the following approximation 
for the line-widths within the WBL approximation, 

* A« m (t,r)«A^ m . (D2) 



At time t = the entire fully connected system (D + L + 
i?) is in its ground state with the chemical potential /x°. 
Afterwards the external potential is switched on, result- 
ing in homogeneous time-dependent level shifts Ae a (t) 
for the lead a (L or R). Hence, for t, t > we have 



S a,nm( T ,*) = h nk a { T ) h Km{t) 9k a i T ^) 

k a Ea 

= hnk a (T) hk a m{t) 

x [i ,r(£fe) e^ (t ~ T) e*£ Ae " {i) Al 



+ 30 



^-rlEtt'Wtf'W, (D4) 



leD 



while for r < and t > 0, the counterparts of (|D3[) and 
(jD4[) arc as follows, 



Sa, nm (r,t) = ^ h nka (T)h kam (t)g< a (T,t) 

k a £a 

= ^ ^«fe Q ( T ) h k a m(t) 

x [< / Q (e£) e* e * ( * _T) e i J? A£ ° (f) dl 
= — o * Jo Ae ° (*) *A" 



+ 00 



/ a (e) e^-^de \ , (D5) 



Gr im (t,r) = ^^7 ) (t)GF m (0,r). (D6) 



leD 



Here the effective propagators for the reduced system, 
are defined as 

U (± \t) = expj ±i J h D {r)dT± Atj, (D7) 



where A = I]a=Lfl AQ - % inserting Eqs. (fDTj) - (|D6|) 
into Eq. (JSj) , the dissipation term Q Q is simplified to be 

QWBL {t) = K - {t ) + {K-,a D {t)}, (D8) 

where the curly bracket on the RHS denotes an anticom- 
mutator, and K a {t) is a Hermitian matrix: 



K a (t) = P a (t) + [P a (t)} 1 . 



(D9) 



Here P a (t) involves an integration over the entire real t- 
axis, which is then decomposed into positive and negative 
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parts, denoted by P& '(t) and P Q \t), respectively. We 
thus have 



/-roc 
dTG r D (t,T)X<(T,t) 
-oc 

= PtH^ + P^it). (D10) 
Pt\t) and Pi +) (t) are evaluated via 

Pt\t) = - f dTG r D (t,T)Z<(T,t) 
J —oo 



2i 

n 



exp<{i/ ke a {T)dT\U ( -\t) 



dee 1 ' 



e - h D {Q) + iA 



A Q , (Dll) 



and 





Pi +) (t) = --f deWtH 

w J-oo 



e,t) 

x / dTWi +) (e,T)A a , (D12) 

respectively, where 



However, the evaluations of Eqs. (|D12|) — (ID13P are found 
extremely time-consuming since at every time t one needs 
to propagate W^(e,t) for every individual e inside the 
lead energy spectrum. It is thus inevitable to have a 
simpler approximate form for Pi + \t) with satisfactory 
accuracy retained. Note that Eq. (|D12j) can be reformu- 
lated as 



- de [ dr 

n J-oo Jo 



x e ~i J*[hn(i)-tA-Ae°'(t)-e]dt ^ (TJ14) 



For cases where steady states can be ultimately reached, 
Ae a (f) and hn(t) become asymptotically constant as 
time t -> +00, i.e., Ae a (t) -> Ae Q (oo) and h D (t) -> 
h D (oo). Therefore, the steady state Pi + \oo) can be 
approximated by substituting Ae"(oo) and ho (00) for 



Ae a (t) and h D (t) in Eq. (|D14j) . respectively. 



2 r» 



de / dr 

n J-oo Jo 
x g— i[fcn(oo)— iA— Ae"(oo)— e](t— t) 

~- f" {' 

I" J-00 1 



-z [ftn (00) — iA— Ae" (00) — e]t 



X 



e - h D (oo) + iA + Ae a (oo) 
It is obvious from Eq. (|D14|) that 



A". (D15) 



PW(Q)=0. (D16) 
Thus Pa (t) for any time t between and +00 can 
be approximately expressed by adiabatically connecting 
Eq. (015]) with ||D16|) as follows, 

p(+)(t) « -— j^I - e -if£lhn(T)-iA-Ae°'(T)-e]dT^ 

de 



-h D {t) + ik + Ae a (t) 



A a . 



(D17) 



Both Eqs. (jDi"4)) and ||DT7|) lead to the correct P Q (oc) 
for steady states, 



P a (oo) 



* r * 

7T 



1 



e - /id(oo) + iA + Ae a (oo) 



A a ,(D18) 



If the external applied voltage assumes a step-like form, 
for instance, AV a (t) = -Ae a (t) = AV a (l - e"*/ a ) with 
a — > + , and hjj(t) is not affected by the fluctuation 
of <7£>(i), Eq. piVp wou ld recover exactly Eq. (|D14|) . 
In other cases, Eq. (|D17j) provides an accurate and ef- 
ficient approximation for Eq. (|D14[) . so long as AV a (t) 
do not vary dramatically in time. Since the integration 
over energy in Eq. (|E)17[) can be performed readily by 
transforming the integrand into a diagonal representa- 
tion, Eq. (|D17|) is evaluated much faster than Eq. (|E)14[1 . 
Due to its efficiency and accuracy, Eq. (|E)17|) is then com- 
bined with Eqs. (|D8[) — (jDlip to calculate the dissipation 
term Qf BL , and thus recovers Eq. ^ of SecHTTBl 
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